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I. INTRODUCTION 



Quarkonium inclusive annihilation decays are historically among the earliest applications 
of perturbative quantum chromodynamics (QCD) [1-4]. At present, it is widely accepted 
that these quarkonium inclusive decay processes can be systematically described by the 
nonrelativistic QCD (NRQCD) factorization approach [5], which is based on the effective- 
field-theory formalism and directly linked with the first principles of QCD. 

In NRQCD factorization approach, the inclusive decay rate can be systematically ex- 
pressed as the sum of product of short-distance coefficients and the NRQCD operator matrix 
elements. The short- distance coefficients encode the hard effects of quark and antiquark an- 
nihilation at length scale of order 1/m (m denotes the mass of heavy quark), which therefore 
can be accessed by perturbation theory owing to asymptotic freedom of QCD. In contrast, 
the NRQCD matrix elements are sensitive to the nonperturbative dynamics that occurs at 
distance of 1/mv or longer (v denotes the typical (anti-)quark velocity inside a quarkonium). 
An important feature of these nonperturbative matrix elements is that they are universal, 
and satisfy a definite power-counting in v. This attractive feature endows the NRQCD 
factorization approach with a controlled predictive power. 

Among the quarkonium annihilation decay processes, the simplest and cleanest are the 
quarkonium electromagnetic decays, exemplified by vector quarkonium decay to a lepton pair 
and pseudoscalar quarkonium decay to two photons. Both of these two processes have been 
comprehensively studied in theory and experiment for decades. In this work, we will address 
the 0(a s v 2 ) correction to the latter process in NRQCD factorization context. We note that 
investigations on this process from other approaches are also available (e.g., see [6-8]). 

We will generically label a pseudoscalar quarkonium by t/q, where Q can stand for the c 
or the b quark. Through relative order- 1> 4 , the NRQCD factorization formula for the decay 
rate of t\q — > 77 reads [9, 10]: 



where ip and y} represent Pauli spinor fields that annihilate the heavy quark Q and heavy 
antiquark Q, respectively. 
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where eeg denotes the electric charge of the heavy quark Q. The 0(a s ) correction to F( 1 S'o) 
was first computed in [11-13]. An incomplete 0(a 2 ) correction to this coefficient has also 
been available about a decade ago [14], which indicates an uncomfortably large negative 
correction. The coefficient G( 1 S , ) is associated with the order-f 2 matrix element, whose 
tree-level value has been known long ago [15]. Recently, the tree-level coefficients H( 1 Sq) 
associated with the order-w 4 matrix elements are also available for the first time [9, 10]. 
Initially only the combination of two 0(v 4 ) coefficients was given [9], as shown in (2c). 
Later Ref. [10] was able to determine these two coefficients separately: iJ 1 ( 1 S'o) = 20ire Q a 



In recent years there have also been many phenomenological investigations on the nonper- 
turbative NRQCD matrix elements for the pseudoscalar quarkonium state [16-18]. Most of 
them focus on the matrix elements at lowest order (LO) and next-to-leading order (NLO) in 
v 2 . Comparing equation (1) truncated at order v 2 with the measured decay rate of i] c — > 77, 
some authors are able to fit the first two NRQCD matrix elements of rj c and find they roughly 
obey the velocity counting rules [16, 17] 1 . 

By far, the 0(a 2 v °) contribution [14] and O(a s v 4 ) contribution [9, 10] to the decay rate 
of tjq — > 77 are known. Assuming a s (m) ~ v 2 , one may naturally wonder what is the actual 
size of the 0(a s v 2 ) correction. In order to answer this question, one needs first know the 
0(a s ) correction to the short- distance coefficient G^So). It is the purpose of this work to 
compute this correction through the perturbative matching method. 

The remainder of this paper is organized as follows. In section II, we outline the per- 
turbative matching strategy which can be utilized to deduce the NRQCD short-distance 

1 The measured partial width of rj c — > 77 seems to have changed significantly in past few years, though 
the full width of r\ c almost remains intact. For example, the branching fraction of this decay channel was 
reported to be (2.4±J;|) x 1CT 4 in PDG 2008 edition [19]. But this value has reduced to (6.3 ±2.9) x 10~ 5 
in the latest PDG 2010 edition [20]. This may cast some shadow on the reliability of the fitted values for 
rj c NRQCD matrix elements by using older data [16, 17]. 
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coefficients for our process. In section III, we elaborate on some technical issues encoun- 
tered in calculating the 0(a s ) correction to QQ^So) — > 77 with the covariant projection 
approach. In particular, we specify the prescription of 75 in dimensional regularization 
adopted in this work. We also mention some technical ambiguities about extracting the 
S'-wave amplitude. In section IV, we then employ the covariant projection technique to 
compute the QQ^So) — > 77 amplitude to NLO in a s . The corresponding 0(a s ) calculation 
in the NRQCD side is presented in section V. In section VI, by comparing the NLO QCD 
amplitude and the respective NRQCD amplitude, we determine the first two short-distance 
coefficients in (1) through order a s . Finally in section VII, we present a brief summary. The 
appendix is devoted to enumerating the analytic expressions of the one-loop scalar integrals 
encountered in the calculation of the NLO QCD correction to QQ^Sq) — > 77. 



II. THE METHOD OF DEDUCING THE SHORT-DISTANCE COEFFICIENTS 

In many effective field theories, it is a standard procedure to determine the short-distance 
coefficient through the matching procedure. We will follow this orthodox method in this 
work. Our calculation is in close analogy with Ref. [21], where the 0(a s v 2 ) correction for 
J/ip — > e + e~ has been deduced also using the matching approach. 

It is worth mentioning that the alternative, even more efficient way of deducing the 
short-distance coefficients exists, exemplified by the method of region developed by Beneke 
and Smirnov [22]. This approach allows one to dissect the loop integral for quarkonium 
annihilation process into four distinct regions: hard, soft, potential and ultrasoft. The short- 
distance coefficient only receives the contribution from the hard region, while the NRQCD 
effective theory characterizes the dynamics from the three low-energy regions. To determine 
the short- distance coefficient, one can either directly compute the hard region contribution 
order by order in v expansion, or by first calculating the full QCD diagram, then subtracting 
the contributions from three low energy regions. The second strategy turns out to have some 
great technical advantage, which enables one to deduce a class of relativistic corrections to 
arbitrary orders in v 2 with ease. Following this spirit, there have recently been progresses 
in inferring the all-order-in-f 2 corrections to J/i/j — > e + e~ [23] and B c — > ez/ e [24]. 
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A. NRQCD factorization formula for the decay amplitude and width 



Let us label the momenta of t]q and two photons by P, k\ and k 2 , and the polarization 
vectors of two photons by e 1 and e 2 , respectively. The Lorentz and parity invariance dictates 
the amplitude of t]q — > 77 uniquely of the structure e a p lJLU P a k^e* 1 ^e* 2 u . For the problem at 
hand, it is most natural to work in the tjq rest frame. In this frame, the amplitude is then 
proportional to the kinematic invariant kx • e* x e 2 . 

When considering quarkonium electromagnetic annihilation decay, one can directly invoke 
NRQCD factorization at the amplitude level. We need only retain those NRQCD color- 
singlet operator matrix elements that connect the vacuum to the r/q state. To the order of 
desired accuracy, one expects that the following factorization formula holds: 

A[ VQ 77] = k ■ el x s* U(0|xVl^) + — 2 (Q\x\-i D) 2 ^q> + 0(v 4 )} , (3) 

where q (i — 0, 2) signify the short-distance coefficients associated with the NRQCD ma- 
trix elements at LO and NLO in v 2 . It is especially convenient in this work to adopt the 
nonrelativistic normalization for t/q state in both sides of (3). In addition, we have explic- 
itly factored out the kinematic invariant, which is represented by a dimensionless Lorentz 
pseudoscalar, 

ki • e{ x e* 2 = --!Le a ^P a kfere* 2 \ (4) 

VQ 

in the right-hand side of (3), where ki = ki/|ki| is a unit vector. The separation of this 
kinematic factor renders q {i — 0, 2) to bear a very simple form. 

Squaring the amplitude in (3), summing over the polarizations of photons, and integrating 
over the phase space and accounting for the indistinguishability of two photons, we can 
express the decay rate of r/g — > 77 as 

r I*^ = k 1 (2% - ^ - U XM» -> 771I 2 

1 c (0|xVl^) + 7 ^(0|x t (-fD)^|r ?Q ) + ---| 2 , (5) 
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where the formula ^ Pol 



ki • e\ x e* 2 



2 



= 2 has been used. 

Comparing (5) with (1), one can express the short-distance coefficients F^Sq) and G r ( 1 S , ) 
that appear in (1), the standard NRQCD factorization formula for the decay rate, in terms 
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of Co and C2- 



o 

Tfi 

FCSo) = ^|c | 2 , (6a) 



m 2 



G(%) = — Re[c c*}. (6b) 

The goal of this work is to determine G( 1 5 , ) through the order a s . Therefore, we first need 
determine both of the coefficients c and c 2 to order a s . 



B. Matching at the amplitude level 

To determine the values of Cq and c 2 , we follow the moral that these short-distance 
coefficients, encapsulating the hard quantum fluctuations that occur at the length scale 
~ 1/m, are insensitive to the long-distance hadronic dynamics. As a convenient calculational 
device, one can replace the physical i]q meson by a free QQ pair of the quantum number 
1 Sq \ so that both the full amplitude, AlQQ^S^) — > 77], and the NRQCD operator matrix 
elements can be accessed by perturbation theory. The short- distance coefficients q then 
can be solved by equating the QCD amplitude A and the corresponding NRQCD amplitude 
A NR q CD , order by order in a s . This procedure is commonly referred to as perturbative 
matching. Analogous to (3), one can write down the pertubative matching formula for the 
QQ^S^) — > 77 process: 

A[QQCS l ^ ] ) 77] = k! • el x e* A NRQCD , (7a) 
Anrqcd = c (0\x^\QQCS [ o ] )} + ^(0| X f (-f D) 2 ^^^ 1 ^ 1 )) + • • • • (7b) 

In Eq. (7), we again adopt the nonrelativistic normalization for the Q and Q states in the 
computations of the full QCD amplitude and the NRQCD matrix elements. 

One can organize the full amplitude A in powers of the relative momentum between Q 
and Q, denoted by q: 



^[QQ^St 1 ) ^ 77] =ki- el xe 



+ ^2 + 0(q 4 ) 



(8) 



For convenience, we have again factored out the kinematic invariant ki • e* x e\ in the 
amplitude. To the desired accuracy, one can truncate the series at 0(q 2 ), with the first 
two Taylor coefficients denoted by £? and srf 2 - To our purpose, we need compute both s^i 
through NLO in a s . This will be conducted in Section IV. 
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The NRQCD amplitude A NR q CD in (7), or equivalently, the NRQCD vacuum-to-QQ 
matrix elements, need also be worked out in perturbation theory through 0(a s ). The 
encountered NRQCD matrix elements at LO in a s are particularly simple 2 : 

(0\x^\QQCS [ o ] )) {0) = VlNc, (9a) 

(o\xK-^) 2 MQQCs l o ] )) (0) = v^q 2 , (9b) 

where the factor y/2N c is due to the spin and color factors of the normalized QQ^Sj^) state. 
The computation of these matrix elements to 0(a s ) will be addressed in Section V. 



III. TECHNIQUES ABOUT COMPUTING FULL QCD AMPLITUDE 

In this section, we outline some necessary techniques about calculating the amplitude for 

gg( 1 4 1] ) 77- 



A. Kinematic setup 

Let p and p represent the momenta carried by Q and Q, respectively. It is customary to 
decompose them in the following form: 

p=\p + q, (10a) 

P=\P-Q, (10b) 



where P is the total momentum of the QQi^S^) pair with invariant mass \[P 2 = 2E, q is 
the relative momentum. Enforcing Q and Q to stay on their mass shells, one requires that 



E = a/ m 2 — q 2 and P ■ q — 0. 

In the rest frame of the QQ pair, which is our default choice, the explicit forms of all the 
momenta are given by 

P» = (2£,0), (11a) 

<f = (0,q), (lib) 

jf = (E,q), (11c) 

j? = (E,-q). (lid) 



2 Unless otherwise stated, throughout this work we use the superscripts (0) and (1) to indicate the LO and 
NLO contributions in a s , and the subscripts and 2 to represent the LO and NLO contributions in v 2 . 
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Hence the total momentum P becomes purely timelike, while the relative momentum q is 
purely spacelike, and one has q 2 = — q 2 and E = ^m 2 + q 2 . In this frame, the momenta 
carried by both photons have a magnitude of E. 
For latter use, we define two velocity variables: 



|qj_ 

~ E ' 



(12a) 



v = M (12b) 



m 

We will distinguish these two variables even in the nonrelativistic limit. 



B. Covariant projection approach 

1. Projection of spin- singlet QQ state 

We start with the quark amplitude Q(p)Q(p) — >■ j(ki,ei) + 7(^2^2), with the momenta 
of Q and Q defined in (10): 

u(p)Tv(p) — Tr \y(j))u(p)T] . (13) 

Here T denotes a matrix in Dirac-color space. 

To proceed, we need first to project the amplitude (13) onto the spin-singlet color-singlet 
Q(p)Q(p) state, by replacing the v{p)u{p) with a suitable projection matrix. The projector 
that is valid to all orders in q for the spin-singlet color-singlet channel, denoted by U^\p,p), 
was first derived in [9]: 



_^ 11 1 

U { i\p,p) = ^2u(p,s 1 )v(p,s 2 )(-,s 1 ',-,s 2 \00)®-j: 



c 

Sl,S2 

' (p' + m)(P+2E) l5 {p--m) <g> -^=, (14) 



8V2E 2 (E + m y ^N~c 
where l c is the unit matrix in the fundamental representation of the color £77(3) group. 
The above spin-singlet projector is derived by assuming the nonrelativistic normalization 
convention for Dirac spinor. Applying this projector to (13), one obtains the amplitude for 
a spin/color-singlet QQ pair annihilation decay into two photons: 

A**[QQ 77] = Tr{u?(p,p)T\, (15) 



where the trace is understood to act on both Dirac and color spaces. 



2. Projection of S-wave amplitude 



In the amplitude in (15), the QQ pair is warranted to be in the spin-singlet state, but 
not necessarily in the S-wave orbital-angular-momentum state. To project out the S-wave 
amplitude, one needs average the amplitude „4 smg over all the directions of the relative 
momentum q in the QQ rest frame. Afterwards one expands the resulting S-wave amplitude 
in powers of q 2 and truncate the series at the desired order. 

A question may concern us immediately- should the relative momentum q be taken as a 4- 
dimensional or a D-dimensional vector upon S-wave angular averaging? Since for our process 
the matching can be conducted solely at amplitude level, there is no a priori criterion to tell 
which treatment is superior. At first glancing, our incapability to pick up a "unique" and 
"correct" scheme may look troublesome, since different treatments may, conceivably, lead to 
different answers for the full QCD amplitude once beyond LO in q 2 expansion. Nevertheless, 
it is possible that both schemes are equally acceptable, provided that ultimately they yield 
identical short-distance coefficients. In this respect, the ambiguity about the dimensionality 
of q may turn into a virtue, in that it can serve as a useful consistency check against our 
calculation. In Section IV, we will separately treat q^ to be 4-dimensional (dubbed g 4 
scheme), and .D- dimensional (labeled qn scheme) 3 . 

Since we are only concerned with the first-order relativistic correction, we will follow a 
standard shortcut to extract the S-wave amplitude rather than literally perform the angular 
integration. In g 4 scheme, we first expand the spin-singlet amplitude A sing in g M through the 
quadratic order, then make the following replacement: 

q 2 

q»q v ->■ — W V {P) (q 4 scheme), (16) 

3 

where 

W\P) = -g^ + —. 

Subsequently we would be able to identify (i — 0, 2) as indicated in (8). 

3 Since q° = in the rest frame of the QQ pair, it seems more natural to view q^ as either a 3-dimensional 
or a D — 1-dimensional vector. Nevertheless, in the lack of confusion, we will stick to the terms #4 and 
q D scheme. 
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If q 11 is assumed to be a D-dimensional Lorentz vector, one can follow the same step as 
described above, except one make the following substitution: 

q>lqU ~* 77ri nM ^ P ) (g D scheme). (17) 

We stress that in both schemes, the momenta P, ki, k 2 and polarization vectors e±, e 2 
are always assumed to reside only in physical spacetime dimensions. 

C 75-prescription in Dimensional Regular izat ion 

When applying the covariant spin-singlet projector as given in Sec. Ill B 1, one needs to 
deal with the trace of 75 with a string of Dirac 7-matrices. This will pose one notorious 
problem, that a definite prescription of 75 must be specified if the spacetime dimension D 
is deformed from four. 

We first point out that, for our process, in principle there exists no any technical subtlety 
about 75 in the q± scheme. One can always choose to first calculate the quark amplitude 
QQ 77 through NLO in a s . As usual, dimensional regularization (DR) can be chosen to 
regularize both ultraviolet (UV) and infrared (IR) divergences. After the loop integration 
is done, one will end up with the T-matrix in (13) that only depends on the external 
kinematic variables, e.g. the momenta of quarks and photons, as well as the polarization 
vectors of photons, which are all 4-dimensional objects. Upon projecting out the spin- 
singlet amplitude, it is obviously legitimate to use the standard 4-dimensional trace formula 
involving 75 in (15). 

In the qu scheme, the 75 problem cannot be circumvented even if one first carries out the 
loop integration for the quark amplitude, because the quark momenta p and p appearing in 
the spin-singlet projector (14) can now penetrate into unphysical dimensions. In this case, 
the rule about D-dimensional trace operation involving 75 in (15) must be specified. 

In practical computation, it is simpler to apply (15) prior to carrying out the loop inte- 
gration. Since the internal fermion propagators and vertices are all D-dimensional objects, 
and all of them will enter into the trace, it becomes compulsory to specify the prescription 
of 75 in DR, for both q D and g 4 schemes. 

In literature, there are two popular prescriptions about 75 in DR, the naive dimensional 
regularization (NDR) [25, 26] and 't Hooft-Veltman dimensional regularization (HVDR) [27, 
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28]. In the former prescription, one assumes {75, 7 M } = for all fj, — 0, 1, • • • , D — 1. In the 
latter, one explicitly constructs 75 = 27°7 1 7 2 7 3 , which anticommutes with 7^ for /i — 0, 1, 2, 3 
but commutes with 7^ for /x = 4, • • • , D — 1. 

In an arbitrary dimension, the definition Tr[7 5 7 M 7 !/ 7 Q 7 /3 ] = — 4ie / ^ ct/3 and {75, 7^} = are 
incompatible [28], therefore 75 in NDR is an ambiguous object. In order to obtain consis- 
tent predictions in this scheme, one must impose some additional rules, e.g., to give up the 
cyclicity property of trace and to place 75 in a fixed position called "reading point" [26] . Nev- 
ertheless, for its technical simplicity, the NDR scheme has been widely utilized in computing 
the NLO QCD corrections to quarkonium decay and production processes [29], though most 
of which are at the LO accuracy in v only 4 . 

In contrast to NDR, the HVDR scheme turns to be a mathematically consistent scheme, 
in which 75 is a well-defined and unique object. For example, the HVDR scheme can 
automatically guarantee to recover the celebrated Adler-Bell-Jackiw chiral anomaly, while 
in the NDR scheme, some ad hoc prescription has to be imposed to achieve this. 

In HVDR, the 7-matrices in D dimension obey the following anticommutation algebra: 

{7^, Y} = 2<T , {75, 7 M } = 2^757*- (18) 

g^ u denotes the projection of the metric tensor onto the unphysical dimensions, which equals 
g^v for /i, v — 4, • • • , D — 1, and equals otherwise. Some useful relations about this tensor 
are gfi = D — 4, g^g" = g^aQv = 9^- A nuisance of the HVDR scheme is that, since the 
first four dimensions are singled out as special, Lorentz covariance has been sacrificed in the 
intermediate stage. Moreover, one may get the impression that the messy anticommutation 
rule for 75 in D dimension would render the practical calculation a formidable task. 

In this work, in favor of its internal consistency, we choose to work with the HVDR 
scheme. It is necessary to spell out the recipe of D-dimensional trace operation involving 
75 in this scheme. In arbitrary spacetime dimension, the trace of a 75 with odd number of 
7-matrices always vanishes. Starting from (18), West has derived a recursive formula for the 

4 It has also been pointed out in [29] , that the 75-prescription implicit in the threshold expansion technique 
developed by Braaten and Chen [30, 31] is essentially equivalent to NDR. 
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trace of 75 with an even number of 7-matrices [32]: 




(19a) 



(19b) 



where the Levi-Civita tensor is a 4-dimensional object in HVDR, i.e., g^^p-y = 0. 

One attractive point of West's trace formula is that, it involves only the 4-dimensional 
antisymmetric tensor together with the D-dimensional metric tensor (not the evanescent 
metric tensor g^l). As a consequence, this recursive algorithm can be readily implemented 
in the Mathematica packages specialized to high energy physics, such as FeynCalc [33]. 

Equation (19) is valid in any dimension, of course also in D = 4, however it is superficially 
much more involved than the familiar 4-dimensional trace formula 5 . 

To make use of West's formula when projecting out the spin-singlet amplitude, we can 
employ the cyclicity of trace to move 75 in (15) to the leftmost, since this property persists 
to be a valid operation in the HVDR scheme: 



and the color trace has been implicit. 

For obvious reason, we would not bother to use the anti-commutation relation of 75 to 
realize this goal. 

IV. QCD AMPLITUDE OF QQ^S^) -> 77 THROUGH NLO IN a s 

In this section we employ the covariant projection technique described in the preceding 
section to compute QQ^Sq^) — > 77 through NLO in a s . The HVDR scheme will be used 
throughout. 

5 For example, in the HVDR scheme, the trace of 75 with six 7-matrices will result in 15 terms, in contrast 
to the 6 terms that one would directly obtain in 4 dimension. 



A si ^(QQ 77) 



1 



Tr<^ 75 (p - m)T (j/> + m)(p+2E) 



(20) 



8V2N~ C E 2 (E + m) 



12 



AAAAAAAA, 



VWXAAAA/ 

a) 



-AAAAAAAA/ 



JVWWWV 



6) 



AAAAAAAA, 



AAAAAAAA, 



AAAAAAAA, 



AAAAAAAA, 



AAAAAAAA, 



AAAAAAAA, 



d) 



AAAAAAAA, 



-AAAAAAAA, 

«0 



FIG. 1: Feynman diagrams for QQ( 1 5j, 1] ; 



77 through 0(a s ). For simplicity, the crossed dia- 



grams have been suppressed. 



A. Tree-level amplitude and matching coefficients 

There are two O(a s ) dia grams for QQ — > 77, one of which is illustrated in Fig. la). The 
corresponding LO T-matrix reads: 

; jzf— + m — p + ^1 + m 



T (o) 



,2 2 



-ie e 



Q 



lr. 



(21) 



-2p ■ k\ —2p ■ k\ 

If q lives only in physical dimensions, one can directly substitute (21) into (15) to project 
out the spin-singlet amplitude, and use the 4-dimensional trace formula to obtain [9]: 

mE 2 



A s[a ^°\QQ -> 77) = e^lVmkt -elxet 



(22) 



' 2 E*-(h- q y 

If q is instead allowed to leak into the unphysical dimensions, one needs substitute (21) 



into (20), and utilize (19) to carry out the D- dimensional trace 
A^(QQ^ 77) =e 2 4\/T 1 



x 



(E + m)(E* - (h ■ q) 2 ) 
E(E + m)e^k la + e^k la q^q u - e ua ^k la q^ + e^q a h ■ q 



Pf3 £ lu £ 2u- (23) 



We have dropped those terms linear in q, which trivially vanish due to the constraints 
P ■ Si = and P° = 2k®. Note this expression is much more complicated than (22), though 
they should be exactly identical when q is a 4-dimensional vector. 

The Levi-Civita tensor, k\, S\ and e<i are all 4-dimensional quantities. For any q vector in 
(23) that contracts with them, only its first four components can contribute, that is, one can 
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make the replacement q a — >■ q a = (g^ — 9p)qp, without affecting the answer. The unphysical 
components of g, g a = g^qp, contribute to (23) only implicitly through the factor E. 

Expand (22) and (23) in q, then apply (16) and (17), one can pick up the corresponding 
tree-level S'-wave amplitudes in the g 4 and qr> schemes. It is then straightforward to identify 
as introduced in (8): 



/(o) 



,(o) 



/(0) 



'2N r 



m 



54 scheme 

= -^2N Q . 

q D scheme (D - l)m 



(24a) 
(24b) 

(24c) 



Not surprisingly, the 0{y 2 ) S'-wave amplitudes differ in g 4 and q® schemes. 

The NRQCD matarix elements at LO in a s have been given in (9). According to (7b), 
one can write down the Born-order perturbative NRQCD amplitude: 



A 



(0) 

NRQCD 



+ 4°V + 



(25) 



Combining (7), (8), (24) and (25), one easily recognizes the O{a s ) short-distance coeffi- 



cients cf ) = ^ } ly/2N c (i = 0, 2): 



,(o) 



.(0) 



.(0) 



c (0) 



4:ireQa 



m 



<74 scheme 



qz> scheme 



8irena 



3m 



3-2e 



m 



(26a) 
(26b) 

(26c) 



where as usual, D = 4 — 2e. The coefficient in q D scheme differs from that in q 4 scheme 
by an (9(e) constant. At first sight, it seems no need to retain this extra piece since no any 
divergence emerges at this order. However, as will become clear later, this 0(e) piece plays a 
key role for ultimately obtaining the scheme-independent 0(a s ) short-distance coefficients. 



B. QCD amplitude at NLO in a s 

We proceed to compute NLO QCD correction to the QQ^Sjp) — > 77 process. At 0(a s ), 
there are eight one-loop diagrams, including two self-energy diagrams, four triangle diagrams 
and two box diagrams. Half of these diagrams have been shown in Fig. 16) through Id). 
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Each individual NLO diagram may contain UV or IR divergences. For example, the 
self-energy and triangle diagrams contain UV divergences, and the box diagrams possess 
IR divergence. We will choose DR as a convenient regulator to regularize both types of 
divergences. Since our concern is to calculate the gauge-invariant on-shell amplitude, for 
simplicity we will work with Feynman gauge in this Section. 

In accordance with the LSZ reduction formula, we need multiply the tree-level amplitude 
in Fig. la) by the residue of the heavy quark propagator at its pole, Zq. This contribution 
is represented by Fig. le). In Feynman gauge, the residue is given by 



J Q 



2 



4tt Ve uv e m m 2 



N 2 — l 

where 7 £ is the Euler constant, and Cf = 2 n c ^ s the Casmir for the fundamental represen- 
tation of the SU(N C ) group. 

In addition, we also need replace the bare quark mass in the quark propagator in Fig. la) 

by 



l-%^(^-3 7 , + 31n^ + 4) + (9(a 2 



(28) 



m w = m , 

An \ euv m 

It is straightforward to write down the respective T^-matrix for each NLO diagram in 
Fig. 1. We then substitute them into (20), and for simplicity, use (19) to carry out the 
.D-dimensional trace prior to performing loop integration 6 . 

Practically, we resort to the Mathematica package FeynCalc [33] to accomplish the 
abovementioned trace operation, because West's formula (19) is its built-in algorithm for 
calculating the trace involving 75. We continue to use FeynCalc to reduce all the en- 
countered one-loop tensor integrals to the one-loop scalar integrals. It turns out that only 
a couple of two-point, three-point scalar integrals and one four-point scalar integral are re- 
quired. For reader's convenience, the closed-form expressions for these scalar integrals have 
been collected in Appendix A. 

By far we have obtained the analytic expression for A smg ^\ the 0(a s ) spin-singlet ampli- 
tude for QQ — > 77. We proceed to expand it to second order in g, then apply (16) and (17) 
to extract the corresponding S'-wave amplitudes defined in (8), in both g 4 and q D schemes, 



Since 75 in the HVDR scheme is mathematically unambiguous, reversing the order of trace and loop 
integration would not affect the final result. 
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TABLE I: The individual contributions to the QQ^S^) ->• 
different classes of diagrams in Fig. 1. A common factor y/2N, 

the Born-order results, while a common factor y/2N, 



27 amplitude, £/q and 3/2, from 



has been suppressed for 



AC F eQaa t 



has also been dropped for the 



0{a s ) results. For brevity, we have used the shorthand i = - — je + hi 



QCD diagrams 


04 scheme 


go scheme 




la) 




1 


lb) 


-J-+21n2-f 


lc) 


1 1 

2 euv 


- 2 In 2 + ^ 


Id) 


7T^ , 17T / 1_ 

4v 4v \ eiR 


+ 2 l n(2 „)) + i I i__l 



le 



\ euv em, y 



la) 
lb) 
lc) 
Id) 

le) 



+ - 21n2 - § 

3 e IR 9 

_5^L 1 1J 2 

6 euv 3 eiR 3 



lluV-§ ln2 + fl 
4luV + fln2-*-I 



£ + ^ (-1^ + 2 In(2«)) £ + ^ (-^ + 2 ln(2t,) + 



+ 



l l 

3£"m 

5 1 

6 e UV 



2 In 2 - i 
_ i_i l 

" 3 em 2 



respectively. The intermediate steps are straightforward but cumbersome, and some special 
care has to be paid when dealing with the box diagrams. Fortunately, almost all these 
manipulations can be handled by computer. 

For completeness and for clarity, we tabulate in Table. I the individual contributions 
to the S-wave amplitudes srf^ and sf!^ from each class of diagrams, in both q 4 and q D 
schemes. For each entry, we sum all the diagrams that belong to the same topology class, 
including the crossed ones. The entry "le)" in Table I refers to the 0(a s ) contributions 
from both wave function and mass renormalization, as specified in (27) and (28). 

As can be seen from Table I, the Coulomb singularity arises in the box diagrams, present 
in both srf^ 1 and srf^ • The emergence of this 7r 2 /t> singularity reflects that the real part 
can receive the contribution from the potential region. For completeness, we also explicitly 
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include the contributions of the imaginary part from the box diagrams. This corresponds to 
a configuration that the two non-adjacent Q and Q quarks in the box can become simulta- 
neously on their mass-shells, so that the exchanged gluon can only be potential mode. Con- 
sequently, only the potential region can contribute to the imaginary part, which is plagued 
with the joint iir/v and IR singularities as well as the term nonanalytic in q. We expect 
that the corresponding 0(a s ) NRQCD calculation in next section will exactly reproduce 
such potential-region contributions. 

Another observation from Table. I is that, both the q 4 and q D schemes differ on the 0(v 2 ) 
contributions in each individual diagram. This pattern is in sharp contrast to the O(v ) case. 

Summing up all the individual contributions listed in Table. I, one finally obtains the 
complete NLO QCD amplitude for QQ^S 1 ^) ->■ 77: 



m , AC F e 2 Q aa s ^ 2 ™ - 2 ' 1 ~ 2 



m 

/ 2 (1) = ^2N C F 

m 



tx 1 - 20 tH i it ( 1 , q 

+ 7e + In — ^ 



Av 4v \ e IR 7r/x z 
L , 4tt/x 2 \ 19 7T 2 4 



(29a) 



5n 2 hi n ( 1 , q 2 , 

+ — h — — ( + 7b + In + C q 



(29b) 



Interestingly, both the g 4 and qu schemes now fully agree on the real part of the complete 
0{v 2 ) amplitude, and only differ slightly on the imaginary part. The scheme dependence is 
encoded in C q , an imaginary constant, which equals in g 4 scheme, and equals ^ in q^ 
scheme. As will be seen in Sec. VI, the scheme- dependence of C q is intimately correlated 
with that of in (26). 

We make some further comments on the complete NLO QCD amplitude in (29). The UV 
divergences have been swept through the wave function and mass renormalization. At LO 
in v 2 , the IR divergences in the real part completely cancel when summing up all diagrams, 
which is warranted by the neutral color charge of the S'-wave QQ pair [5]. Nevertheless, 
this cancelation fails to hold at NLO in v 2 , and the presence of a net infrared divergence in 
the real part of (29b) signals that the simple "color-transparency" picture no longer applies, 
and the soft gluons can resolve the geometric details of the QQ pair and may even strongly 
interact with its color dipole. It is easy to find that this IR divergence arises from the 
soft region of the loop diagrams. We anticipate that the corresponding 0(a s ) NRQCD 
calculation in next Section will exactly reproduce this infrared divergence at relative order 



v 2 . 
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V. NRQCD AMPLITUDE AT NLO IN a s 



To match the accuracy of the QCD calculation for QQ(}Sq ) — > 77 in Section IV B, it is 
necessary to compute the perturbative NRQCD amplitude A NR qcd to NLO in a s . 

In (9), the two involved NRQCD matrix elements have been given at LO in a s . In the 
following we proceed to compute their 0(a s ) corrections. For this purpose, the knowledge 
of the NRQCD lagrangian in the heavy quark bilinear sector is required [5]: 



The replacement in the last line implies that the corresponding heavy anti-quark bilinear 
sector can be obtained through the charge conjugation transformation. £ii g ht represents the 
lagrangian for the light quarks and gluons. 

The short-distance coefficients cf, cd, cs, in (30) are associated with the Fermi, Darwin, 
and spin-obit interactions, respectively. To our purpose, suffices it to know their tree-level 
values cf = c D = c s = 1 + 0(a s ). Moreover, in (30) we have neglected all other higher- 
dimensional operators, since the relativistic corrections generated by them are beyond the 
intended 0(v 2 ) accuracy. 

Through the relative order v 2 , the 0(a s ) corrections to the matrix elements stem from 
the one-loop NRQCD diagrams as illustrated in Fig. 2. The first six diagrams 7 represent 
the one- loop correction to the matrix element {Olx^lQQi 1 S^)) , while the last one in Fig. 2 
constitutes the one-loop correction to (0|x^(— | ~D) 2 ip\QQ( l S^)) . 

In passing, we note that our one-loop correction calculation for the vacuum-to-QQi 1 S ) 
NRQCD matrix elements is quite similar to the existing one for the vacuum-to-QQ( 3 Si) 
matrix elements [21]. Since heavy quark spin symmetry is violated at relative order v 2 , we 
cannot completely transplant their results here. Ref. [21] uses Feynman gauge in calculating 
NRQCD diagrams. In contrast, in this section we will work in the Coulomb gauge, since 
it has the virtue that the instantaneous (potential) nature of the temporal gluon becomes 

7 We have not drawn the diagram in which a temporal gluon attaches to a spin-orbit interaction vertex and 
a Coulomb vertex, because its contribution simply vanishes due to spin conservation. 



£ N RQCD = V [iDo + 7^ J i> + i> g^3^ + 2^ ° ' # sB ^ 

+ -^^(D • £ S E - g s E ■ D)V + ^^<r ■ (D x £ S E - g s E x D)^ 

Sm z 8m z 

+ (V> -)• io 2 X \ -)■ -Al) + £ light . 





(30) 
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e) /) 9) 

FIG. 2: The one-loop NRQCD diagrams that initiate the 0(a s ) corrections to the vacuum-to- 
QQi}S^) matrix elements. The big solid circle represents the operator ~x}ip, while the solid square 
represents the operator fD) 2 ^- The dashed line denotes the temporal gluon, and the curly- 
line denotes the transverse gluon. The unlabeled vertex attached to a temporal gluon stands for 
the Coulomb interaction ^ig s T a . For simplicity, we have suppressed the conjugate diagrams of b) 
and c), in which the same operators are inserted on the antiquark line. 

manifest. Of course, insofar as the on-shell matrix element is concerned, the choice of gauge 
is merely a matter of convenience. 

When computing the loop diagrams in Fig. 2, we first integrate over the temporal compo- 
nent k° using contour integration, where k signifies the loop momentum carried by the gluon 
propagator. We then perform the remaining integration over k in D — 1 spatial dimension. 
When the exchanged gluon is transverse, as exemplified by Fig. 2c?), the situation is some- 
what more complicated. Upon integrating over k°, one would end up with two terms, one 
with the residue taken at the pole of the quark propagator, and the other with the residue 
taken at the pole of the gluon propagator. The denominators of each resulting term are in 
general inhomogeneous in powers of v. This problem was initially overcome in a heuristic 
way by expanding the integrand [5]. This recipe was later clarified and systemized in the 
method of region formalism [22] (see also [21]). In this context, the two terms after con- 
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TABLE II: The analytic expressions of each 0(a s ) NRQCD diagram in Fig. 2 (in Coulomb gauge) 
For simplicity, we have stripped off an overall factor of V 2N C c ' F ^ Xs . 



NRQCD diagrams 



Expressions 



2a) 
26) 
2c) 
2d) 
2e) 



+ 7^ + ln^-l) 



- 2 - -0 

uv eir / 



3 V £ uv 



2/) 


3 \ euv e lR / 


2a) 


m 2 l 


7r2 + i7r (-^ + ^ + ln ^) 





tour integration can be identified with the potential and soft region, respectively. For each 
region, loop momentum has a definite power in v so that one can readily homogenize the 
integrand by expanding the denominator accordingly. After this is done, it then becomes a 
straightforward exercise to carry out the spatial integrations in DR. 

According to the LSZ reduction formula, we need multiply the LO matrix element 
(0|xVlQQ( 1 Sq))^ by the residue of the heavy quark propagator in NRQCD near its pole, 
Znrqcd- The 0(a s ) piece of this contribution is represented by Fig. 2f). In Coulomb gauge, 
the residue in DR reads 8 

Znrqcd(?) = 1 - (—--)■ (31) 

<jtt \ e uv em/ 

It is worth noting that, the simultaneous occurrences of UV pole and IR pole at 0(v 2 ) can 
be traced to a logarithmically-divergent scaleless integral in the soft region. 

For clarity, we list in Table II the separate contribution from each diagram in Fig. 2. In 
Fig. 2a), we follow [21] to include an extra source of relativistic correction that arises form 

2 4 

the nonrelativistic expansion of the kinematic energy of the external (anti) quark |L — g^-. 
For the entries corresponding to Fig. 26) and c), we have also included the contributions 
from their conjugate diagrams. It is clear to see that, all the diagrams receive contributions 
from the potential region (odd powers of v). In addition, Fig. 2c?) and /) also receive 

8 The Coulomb-gauge residue has been given in Appendix B of [5], which depends on a hard momentum 
cutoff there. We have converted that cutoff-scheme expression to its counterpart in DR. 
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the nonvanishing contributions from the soft region (even powers of v), which result in a 
logarithmically divergent scaleless integral. Only p • A interaction can result in such a 
contribution. 

Summing all the individual contributions in Table II, we obtain the 0(a s ) corrections to 
the perturbative NRQCD matrix elements: 



«Q(- 5 ?>)>« = (l + i/) + in (-i- + - IE + in 



2v 2 



3 \ e uv em, 



(32a) 



lxt(-lB)^|QQ(^)>a) = V^^ 2 Q) 



7T + 27T h 7s + In 



e m 7r/i 2 



(32b) 



As is evident in (32a), the order-o; s correction to the matrix element {0\x^\QQ( 1 So)} is 
logarithmically UV divergent, which can be traced back to the p ■ A interactions in Fig. 2d) 
and /). Employing the MS scheme to subtract the UV divergence, we then obtain the 
renormalized matrix element 

+ -1E + In4n) 1 (33) 



3 \ e iR 

According to (7b), accurate up to the order a s v 2 , we can express the perturbative NRQCD 
amplitude as 

Anrqcd = (4 0) + c^)(0\ X ^\QQCS [ o ] )) {0) + c^(0\ X ^\QQCS [ o ] ))g (34) 
Jo) , J 1 ) ^ Jo) 

+^^(o|x t (-|D)>|gg( 1 ^))(°) + ^xK-^MQWs^ + • • • , 

with all the involved NRQCD matrix elements available now. It is anticipated that (34) will 
faithfully reproduce the infrared behavior of (29). 



VI. MATCHING THE SHORT-DISTANCE COEFFICIENTS AT NLO IN a 



In section IV A we have already conducted the perturbative matching at tree level, and 
determined the tree-level coefficients cf^ (i = 0, 2) in (26). Substituting their explicit values 
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into (34), and subtracting the contribution of the tree-level NRQCD amplitude ^nrqcd> one 
finds that the NRQCD amplitude at 0(a s ) reads: 



7r 2 + in I — — + 7 S + In 



2 



eiR VTyU 2 



+ C/ • (35) 



Here C g is a scheme-dependent imaginary constant. The origin of this constant can be 
easily traced, which arises from the last term in (34), that is, the 0(e) piece in the 
tree-level coefficient multiplied by the imaginary infrared pole in the matrix element 

(0\xK-^mQQCS [ o ] )} {1) in (32b). 

Remarkably, C q turns out to exactly equal the scheme-dependent constant C q that ap- 
pears in the NLO QCD amplitude st!^ in (29b), which vanishes in q 4 scheme, but equals 
in qo scheme. Clearly, the equality of these two constants is by no means a coincidence. 

Following (7) and (8), one readily recognizes the matching condition at order a s : 

^ ( n^ 2 ( V = A« QCD . (36) 

Plugging the analytic expressions of srf^ (i = 0,2), which are given in (29), and the 
expression of A^q CD , as given in (35), into this matching equation, one observes that the 
infrared and Coulomb singularities, together with the terms non-analytic in |q|, exactly 
cancel from the both sides of (36). This is of course as expected, which is warranted by the 
general construction of the NRQCD effective theory. 

It is a nontrivial check for the consistency of our calculation, that the imaginary parts, 
which receive the contribution solely from the potential region, indeed fully cancel in both 
sides of (36). In particular, it is remarkable that the scheme-dependent imaginary pieces 
also cancel in both q 4 and qn schemes, thanks to the equality C q = C q . This cancelation is 
crucial for one to obtain the scheme-independent short-distance coefficients 9 . 



9 We have made a further test of the scheme-independence of the coefficient In the QCD-side calcula- 



tion, we attempt to use the alternative version of the spin/color-singlet projector [34]: 



which is equivalent to (14) only in D — 4. Following the same procedure as described in the text to 
extract the 1 5q 1 ' amplitude, we obtain identical results as in (29) except for the C q term. As expected, 
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It is then straightforward to solve the order-o; s short-distance coefficients: 

A-Ke 2 Q aC F a s ( n 2 5" 



c (1) - 



c (1) 

c 2 



m 7r 

.2 



87re Q o; CfCt^ / // 



19 3tt 2 
4m 2 + 12 ~ "16" 



(37a) 
(37b) 



3m 7r 

where /x is now interpreted as the factorization scale in the MS scheme. Note these coeffi- 
cients are purely real. 

With the knowledge of the c*p in (37), we can employ (6) to determine the short-distance 
coefficients F and G that appear in (1) through order a s : 



2 r 



JO) 1 2 



= 27relo; 2 



1 + 



+ 2Re(4 0) 4 1) *)+O(a s 2 ) 
C F a 



TT 



(38a) 



GCSo) 



m 



Air' 



-Re 



(0)J0)* . Jl) JO)* 











(0) Jl)* 
c 2 



+ 0(a^ 



8-TTeQa 2 



1 + 



Ci?a s / 57r 



12 4m 2 



(38b) 



TT V 16 

The familiar result for F^Sq) through 0(a s ) is then reproduced. The order-o; s correction 
to G( 1 S'o) is new, which constitutes the main result of this work. 

To visualize the relative importance of these 0(a s ) corrections, one may choose \i = m 
in (38): 

ai s (m) 



FCSo) = 2ire A Q a 
GCSo) = 



1-3.38 x 



STreqa 2 



1 + 0.52 x 



TT 

a 8 (m) 



+ 0(a£) 



TT 



+ 0{a 2 s ) 



(39a) 
(39b) 



The order-o; s correction to F( 1 Sq) has modest effect for r\ c decay to 77. If one takes 
®s{ m c) — 0.3, including this correction will reduce the tree-level value of F by about 30%. 
Nevertheless, the order-o; s correction to G^Sq) seems to have a much minor impact, which 
only enhances the tree-level value of G roughly by 5% for \i = m. Note even the relative 
sign of this correction is uncertain. If one instead chooses /i = 2m, including this 0(a s ) 
correction would reduce the tree-level value of G by about 10%. It seems fair to state that, 
at current level of experimental accuracy, the 0(a s ) correction to G( 1 S'o) has very little 
phenomeno logical impact. 



one again gets C q — C q = in q± scheme as before. However, for the qu scheme in this case, one finds 



that C g = C q = 



17T 
' 1SV ' 
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VII. SUMMARY 



NRQCD factorization approach provides a model-independent framework for calculating 
the quarkonium annihilation decay rates. The nonperturbative NRQCD matrix elements 
are universal, which can be estimated from potential models or lattice QCD simulations, 
or directly fitted from experiments. On the contrary, being process-dependent, the short- 
distance coefficients appearing in the NRQCD factorization formula can be systematically 
improved in perturbation theory. To make a more precise prediction, it is desirable to know 
the short- distance coefficients at higher accuracy in a s expansion. 

In this work, we have computed the 0(a s ) corrections to the short-distance coefficients 
relevant to the process of pseudoscalar quarkonium decay to two photons. Specifically, we 
have deduced the 0(ct s ) piece of the coefficient G^Sq), which is associated with the relative 
order-w 2 NRQCD matrix element in (1). This coefficient is determined through equating 
the full QCD amplitude and the NRQCD amplitude for the quark process QQ^Sq) — >■ 77 
through NLO in a s and v 2 . As a consistency check, it is found that the IR and Coulomb 
singularities, as well as those terms nonanalytic in |q|, are exactly canceled upon matching 
the calculations on both sides. We have presented detailed descriptions for some subtle 
technical issues encountered in the calculation, such as consistent prescription for 75 in di- 
mensional regularization, and the ambiguity about the dimensionality of relative momentum 
when projecting out the S*-wave amplitude. We have explicitly verified that different artifi- 
cial schemes in extracting the S-w&ve amplitude actually lead to the identical short-distance 
coefficients. 

Although the 0(a s v 2 ) correction appears to be phenomenologically insignificant for pseu- 
doscalar quarkonium decay to two photons, it is of interest to investigate the analogous 
correction for the process of pseudoscalar quarkonium inclusive decay to light hadrons. 
Note added in the proof. Shortly after this paper was sumbitted, a similar work has 
appeared in arXiv, which also addresses the 0(a s v 2 ) correction to pseudoscalar quarkonium 
decay to two photons [37]. These authors performed the matching calculation directly at 
the decay rate level, and their final results agree with our (38). 
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Appendix A: Useful one-loop scalar integrals 

In this Appendix we tabulate those one-loop scalar integrals that are encountered in the 
calculation of the NLO QCD correction to the process QQ^S^) — > 77. First we adopt the 
following definitions of the 2-, 3-, 4-point Passarino-Veltman scalar integrals [35]: 

Botrf; ml ml) S J^r(l - e) J d °l _ ' + _ (Ala) 



u 2e 

/ 2 2 2 2 2 2\ r 1 -p/i 



X I^ D ^ (I 2 - m\) /n 1 - ^ , „ \9 „„2\> (Alb) 



1 

»?)((/ + fc) 2 -ml)((l + q 2 ) 2 - mlY 



Do(pi,pipi,pi]S 12 ,S23]mi,m 2 2 ,m 2 3 ,ml) = —^T(l-e) 

X I^ 1 (I 2 - m 2 M + <?i) 2 - ™ 2 W + <? 2 ) 2 - m§)((Z + g 3 ) 2 - ™ 2 ) ' (Alc) 
where the spacetime dimension D = 4 — 2e, g n = ^™ =1 Pi and go = and s^- = (p» + Pj) 2 - 
The ie prescription in each propagator has been tacitly assumed. 
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The desired results are: 

B (-E 2 ± 2h ■ q + g 2 ; 0, m 2 ) = J- + In 4 + 2 - ^|^ i ^=^=,(A2a) 

e uv m J 2(£/ 2 + k\ - q) —m l 

Li 2 f - g X ±2fc1 ' 9 )-^ 
C (0, m 2 , -E 2 + q 2 ± 2h ■ g; m 2 , m 2 , 0) = 2( ^ fci . g) > ( A2b ) 

C (4E 2 , 0, 0; m 2 , m 2 , m 2 ) = ^ [ \\n 2 i±| - £ - i vr In i±| j , (A2c) 
( U4E 2 ,m 2 ,m 2 ;m 2 ,m 2 ,0) = -^-<' UL + In^l (~ In + i ^ - vr 2 



4£ 2 /3 1 VdR m 2 / \ 1-/3 
A)(™ 2 , 0, 0, m 2 ; -£ 2 + g 2 ± 2ifei • g, 4£ 2 ; 0, m 2 , m 2 , m 2 ) = ' 



8E 2 (E 2 Th-q)f3 

— + ln^-21n ^^-ln4^ 2 « 

1 — p J VeiR m z m z 



+2 



Li? I —| - Li 



(A2e) 



1 + 13 J V i + Z 3 , 

All the involved kinematic factors, such as k±, g, E, /3, have already been defined in Sec- 
tion III A. The dilogarithm is defined through 

Ii 2 (x) = - / dt ln(1 ~ f) ■ (A3) 

All the infrared divergent one-loop scalar integrals (up to 4-point) have been classified in 
Ref. [35], which also compiles the corresponding analytic expressions. We thus can directly 
deduce the analytic expressions (A2d) and (A2e) from [35] by making appropriate substi- 
tutions. The analytic expression for the Coulomb- divergent three-point integral, (A2d), is 
well known and can be found in many places, e.g., Eq. (4.14) in [35], or Eq. (B10) in [23]. 
However, if one starts from Eq. (4.47) of [35] to deduce the Coulomb-divergent D function, 
one would unfortunately end up with an erroneous expression for the imaginary part. In 
(A2e), we present the correct analytic expression for the imaginary part, which is derived by 
employing the D-dimensional Cutkosky's rule. Its correctness has been numerically verified 
with the help of the Mathematica package LoopTools [36]. 
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